subroutine set_ic( dfu, uxMat, uyMat, f0vec, fvc0T, it)    ! outputs

  use dimensions

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : kz, vPhase, ux0, uy0, &                      ! real,     phys
                         nz, iWaveFrame, iModx, iMody, iModz, &       ! integer,  phys
                         ModAmp,  &                                   ! real,     phys
                         restart, &                                   ! integer,  num
                         pi, iunit, &                                 ! constants
                         runname

  implicit none

!-- interface input/output --

  real,    intent(out) :: dfu(jv, jz, jy, jx, jDirection)
  real,    intent(out) :: uxMat(jz, jy, jx, jDirection)
  real,    intent(out) :: uyMat(jz, jy, jx, jDirection)
  real,    intent(out) :: f0vec(jv)
  complex, intent(out) :: fvc0T(jv)
  integer, intent(out) :: it


!-- local variables ---

  integer              :: ix, iy, iz, iv, iDirection
  real                 :: MoDinZ(jz), MoDinX(jx), MoDinY(jy) 
  real                 :: v, f0


  f0(v)=( 1/sqrt(2*pi) )*exp(-.5*(v+vPhase*iWaveFrame)*(v+vPhase*iWaveFrame))


!-- create vector "f0vec" --


  do iv=1,jv
     v = vvec(iv)
     f0vec(iv) = f0(v)   
  end do
                 
  fvc0T(:)=f0vec(:)      
  call fft_forward(jv,fvc0T,fvc0T)   !  fvc=FT(dfu)
  fvc0T=fvc0T/float(jv)             !  normalized velocity  FT of equilibriium


!-- read restart or create new --

  if (restart.ne.0) then         ! read restart info

       write(*,*)
       write(*,*) 'reading restart'

       call read_bin(dfu, uxMat, uyMat, it)        

  else                   ! start fresh

       write(*,*)
       write(*,*) 'fresh start'

       it = 0
       dfu=0.0

       do iDirection = 1, jDirection
          uxMat(:,:,:,iDirection) = ux0(iDirection)
          uyMat(:,:,:,iDirection) = uy0(iDirection)
       enddo
 
       if ( iModx ) then
          MoDinX(:) = cos( dkx * xvec(:) )
       else
          MoDinX(:) = 0
       end if

       if ( iMody ) then
          MoDinY(:) = cos( dky * yvec(:) )
       else
          MoDinY(:) = 0
       end if

       if ( iModz ) then         
          MoDinZ(:) = cos(  kz * zvec(:) )
       else
          MoDinZ(:) = 0
       end if


       do ix=1,jx                                                      
          do iy=1,jy
             do iz=1,jz
                do iv=1,jv
                   dfu(iv,iz,iy,ix,:) =  &
                        ModAmp * f0vec(iv) * ( MoDinZ(iz) + MoDinX(ix) + MoDinY(iy))
                enddo
             enddo
          enddo
       enddo

       write(*,*) 'initial data created'
 

  end if  !  end restart test


end subroutine set_ic
